Microclimatic spatial variation buffers local warming in alpine vegetation
B Jiménez-Alfaro, E Fernández-Pascual, A Bueno, C Marcenó….¿?
Abstract (ME&B)
Alpine vegetation is considered to respond to climate change through species distribution shifts with effects on community richness and composition. However, the complex topography of alpine landscapes creates a mosaic of microclimatic niches which might buffer macroclimatic variation, preventing local extinctions of alpine species. The magnitude of microclimatic buffering is context-dependent and may depend on regional factors such as biogeographical settings, topography, or species pools. A key question for understanding the magnitude of this buffering is therefore to test whether spatial topographic variation reflects temporal changes within and between years in different regions. We addressed this question in a long-term vegetation monitoring site in Picos de Europa National Park, a transitional mountain massif between the temperate and Mediterranean macroclimates in northern Spain. In 2008, we established four stations along an alpine landscape for measuring temporal changes in soil temperatures and related changes in plant communities based on two 1 m2 plots replicated around the central point. In 2018, we also measured spatial variation in soil temperatures within each station using 20 additional plots placed at 10 m intervals along with cardinal directions from the central point. We found that spatial variation was higher than 10-year temporal variation for most of the bioclimatic indices we used, but there were exceptions in stations that accumulate less snow, and in less topographically diverse stations. Vegetation variation (Sørensen index) was higher across topographical gradients than across time. These results suggest that microclimatic refugia can compensate for temporal changes, but this compensation might not be homogeneous throughout the alpine landscape. The response of the study system to contemporary climate warming seems more likely to produce a slow re-accommodation of species relative abundances along with topographical variation, rather than local extinctions.
Introduction
xxx
Methods
We conducted all analysis with R (R Core Team 2021), and the code for analysis and creation of the figures and manuscript is available at GitHub (see Data Availability Statement).
Study system
The study was conducted in the central calcareous massif of the Picos de Europa National Park, in northern Spain (Figure 1A). The study area is a biodiversity hotspot for cold-adapted plants in the Iberian Peninsula and a biogeographical hub for Alpine and Mediterranean lineages in Western Europe (XXX). The central calcareous massif occupies c. 50 km2 and supports a high diversity of ecosystems, with alpine vegetation mostly occurring between 1900 and 2400 m a.s.l., with a local species pool of XX species (Jiménez-Alfaro et al. XXX). In 2008, we established a long-term ecological research program for monitoring soil climate and vegetation change. We selected four study sites along a north-south gradient (Figure 1B), reflecting the main macroclimatic variation from Atlantic to Mediterranean influence. The sites were selected to… general design and description of sites… (Borja)
Vegetation and microclimatic data
Temporal survey. In each site, we buried a temperature logger (M-Log5W, GeoPrecision, Ettlingen, Germany; accuracy: +/- 0.1 ºC at 0 ºC, resolution: 0.01 ºC, records each hour) at 5 cm depth in a relatively flat and homogeneous vegetation patch. We surveyed the plant community in two replicated plots of 1 m2 separated 1 m from the logger, identifying species composition of vascular plants and estimating relative cover in %. In each plot, we installed a grid template of 100 microplots (10 x 10 cm each) to sample species frequency according to the standard methodology of GLORIA (XXX). The loggers were replaced by new ones, when needed, to obtain a continuous temperature record from 2008 to 2018. In 2018, we resampled the same plots in the same way to detect potential changes in species presence and frequency. The vegetation data from these surveys, together with the soil temperature collected in the four study sites during 10 years, represent the “temporal survey.”
En este punto podríamos comentar que el long-term monitoring incluye otros cuatro sitios, algo más diferentes, que por ser más diferentes no los consideramos en este estudio, pero que añadimos los datos en supplementario. Así, quedan publicados.
Spatial survey. In 2018, we visited the same areas to study the spatial variation of vegetation and microclimate around the previously sampled areas. Using the long-term temperature logger as the central point, we additionally placed 20 iButtons (Thermochron, iButton, Newbury, UK; accuracy: +/- 0.5 ºC from -10 ºC to +65 ºC, resolution: 0.5 ºC, records each 4 hours) in 20 plots of 1 m2 separated 10 m from each other in the four cardinal directions (Figure 1C). At each of the 20 plots, we identified vascular plants and estimated their relative cover in %. In 2019 we came back to download the data of the iButtons. These data, together with the vegetation data of the iButton plots, is the “spatial survey.”
Bioclimatic indices
We used the microclimatic data of the temporal and spatial surveys to calculate several bioclimatic indices. We did this for each year of the temporal survey (4 sites * 10 years) and each of the 80 spatial plots (4 sites x 20 plots). The spatial data series did not cover a full year, and were missing part of August and September. Thus, for comparison, we also removed this period from the temporal data before calculating the indices. Additionally, for the temporal survey (hourly records), we kept only the same recording hours as in the spatial survey (records each 4 hours).
The bioclimatic indices were the following: (1) bio1 = annual mean temperature; (2) bio2 = mean diurnal range,i.e. the mean of the monthly differences between maximum and minimum temperatures; (3) bio7 = temperature annual range; i.e. the difference between the maximum temperature of the warmest month and the minimum temperature of the coldest month; (4) snow = the number of days of snow cover, considered to be those days in which the maximum temperature was below 0.5 ºC and the minimum temperature was above -0.5 ºC; (5) GDD = growing degrees-day, i.e. the sum of daily mean temperatures for days in which the mean temperature was above 1 ºC; (6) FDD = freezing degrees-day, i.e. the sum of daily mean temperatures for days in which the mean temperature was below 0 ºC. The indices bio1, bio2 and bio7 follow the definitions of WorldClim (Fick and Hijmans 2017). For GDD, we used the 1 ºC threshold for the growing season, following Bürli et al. (2021). For FDD, we transformed the values from negative to positive, so the interpretation would be easier (i.e. higher values equal more freezing). The bioclimatic indices are provided as supplementary material S1.
Before proceeding with further analysis, we conducted a principal component analysis (PCA) (Lê et al. 2008) of the bioclimatic indices for the spatial survey, to identify the main patterns in microclimatic variability. Results of the PCA are shown in supplementary material S2. In summary, the first axis of variation represented a gradient of increasing thermicity, with higher values of GDD, bio1, bio2 and bio7. The second axis represented a gradient of increasing freezing, with higher values of FDD. Snow length was mostly associated to the third axis. For ease of interpretation, we kept only GDD and FDD for further analysis. GDD and FDD were not correlated (r = -0.03) while GDD was correlated to bio1 (r = 0.97), bio2 (r = 0.76) and bio7 (r = 0.45).
Simulation of species extinctions
We modeled species occurrence at the micro-site level as a function of the bioclimatic indices GDD and FDD, using the data of the spatial survey. First, we used non-metric multidimensional scaling (NMDS) with environmental fitting (Oksanen et al. 2019) to visualize the patterns of change in vegetation, and their relation to GDD and FDD.
Then, we used Generalized Linear Models (GLMs, binomial family) to model species presence (1 or 0) in each of the spatial plots (n = 78, we removed two plots with no vascular plants) as a response to the plot’s GDD and FDD. We only did this for species with at least 10 presences in the plots. We only kept models in which at least one of the two bioclimatic indices had a significant effect size (p < 0.05) and for which the value of McFadden’s pseudo R2 (McFadden 1974) was higher than 0.15. For readers unfamiliar with McFadden’s pseudo R2, we point out that it tends to have lower values than R2 in ordinary least squares regression, and that values between 0.2 and 0.4 represent excellent fit (McFadden 1979).
To construct scenarios of climatic changes, we used the temporal survey data (10 years, 2009-2018) to describe situations in which the extreme values of today would become the new average. To do so, we calculated, for GDD and FDD, the maximum and minimum values recorded in the entire period. With the resulting four values, we created four scenarios: hot and snowy (max GDD, min FDD), hot and frozen (max GDD, max FDD), cold and snowy (min GDD, min FDD) and cold and frozen (min GDD, max FDD). We used the GLM models to predict the probably of presence for each species and scenario, considering that a probability of 0 in a given scenario would mean the extinction of the species in said scenario.
Results
Microclimatic variation
Along the 10 years of soil temperature monitoring, we found a general trend of increasing temperature, especially in the southwest and warmest site (Hoyo sin Tierra) (Figure 3). Two of the sites (Los Cazadores and Los Boches) showed a consistent pattern of continuous snow cover during winter (reflected by temperature records around 0º C). In contrast, the two other sites (Hoy Sin Tierri and Hoyo Sin Tierra) showed frost temperatures during most winters (although they may not be constant and intermixed with short snow periods). Maximum temperatures were… Soil temperature corresponding to the spatial surveys also showed high variation within sites (Figure 3). In the four sites we found plots with contrasting patterns of snow cover and frost during winter, reflected by… When looking at the whole microclimatic variation among plots (Figure 4), we found that the four sites were represented by a wide range of temperatures, with the exception of one site (Los Cazadores). The two main axes of variation were mainly related to FDD and GDD, the latter also correlated with…. (Eduardo)
Plant diversity and microclimate
Across the whole study system (temporal and spatial surveys) we recorded 86 species of vascular plants (considering Helianthemum apenninum subsp. urrielense and Helianthemum apenninum subsp. cantabricum as separate species), representing % of the total species pool of the study region. The five most frequent species were Thymus praecox subsp. ligusticus (83 occurences), Anthyllis vulneraria (73), Koeleria vallesiana (59), Minuartia verna (55) and Helianthemum canum (52). Average species richness per 1m2 plot was 13, with the richest plot having 25 species and the poorest one just two.
In the temporal survey (2 visits x 2 plots per site, n = 16) we recorded 42 species in 2009 and 47 in 2019. Of the species recorded in 2009, we did not find again in 2019 the following three: Festuca burnatii, Galium pyrenaicum and Iberis carnosa. Conversely, in 2019 we recorded eight species that we had not seen in 2009: Arenaria purpurascens, Lotus corniculatus, Potentilla crantzii, Sedum album, Sedum brevifolium, Seseli montanum, Silene ciliata and Solidago virgaurea. Of these species that appeared or disappeared from the temporal survey plots, Arenaria purpurascens appeared in 12 10x10 cm cells in 2019, while the rest of the cases amounted to less than 10 cells. The five species with the highest decrease in frequency from 2009 to 2019 (ignoring annual species and species that occurred in less than 10 10x10 cm cells in 2009) were Armeria cantabrica (85% decrease in frequency, present in 13 cells in 2009), Poa alpina (-83%, 18 cells), Salix breviserrata (-48%, 25 cells), Jurinea humilis (-26%, 23 cells) and Ranunculus parnassiifolius subsp. favargeri (-18%, 72 cells). The five species with the highest increases (again, ignoring annual species and species that occurred in less than 10 10x10 cm cells in 2009) were Minuartia verna (+278%, 19 cells), Helianthemum apenninum subsp. urrielense (+87%, 63 cells), Arenaria moehringioides (+85%, 13 cells), Saxifraga conifera (+83%, 24 cells) and Silene acaulis (+39%, 18 cells).
Changes in species frequency were/were not correlated with XX and XX…. (Table 1) (Eduardo). In the spatial survey (20 plots per site, n = 80; including one plot placed in rocks with no vascular plants) we recorded 81 species, 36 of which were absent from the temporal survey plots.
Changes in species composition (Fig. 4) were associated to GDD (NMDS axis 1) to represent a main gradient from the warmest to the coldest microsites, and to FDD (NMDS2) in a gradient of snow versus frost winters. In the compositional space, the sites were poorly differentiated across the four sites (PERMANOVA, adonis) mainly by the effect of one site (Los Boches) with most plots places in the cold and frost margins. Etc…. (Borja).
redo environmental fit
Results of the environmental fit:
print(ef1)##
## ***VECTORS
##
## NMDS1 NMDS2 r2 Pr(>r)
## bio1 -0.99853 -0.05411 0.7391 0.001 ***
## bio2 -0.87326 -0.48725 0.5242 0.001 ***
## bio7 -0.71719 -0.69687 0.1814 0.003 **
## Snw 0.33752 0.94132 0.0606 0.080 .
## FDD 0.27395 -0.96174 0.3295 0.001 ***
## GDD -0.96954 -0.24494 0.7780 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999Correlations between bioclimatic variables:
print(biocor)## bio1 bio2 bio7 Snw FDD GDD
## bio1 1.0000000 0.7632716 0.37483777 -0.17089008 -0.27292081 0.96774811
## bio2 0.7632716 1.0000000 0.79299098 -0.13715459 0.16136989 0.82331833
## bio7 0.3748378 0.7929910 1.00000000 -0.03074698 0.47129871 0.46378408
## Snw -0.1708901 -0.1371546 -0.03074698 1.00000000 -0.30045520 -0.19039540
## FDD -0.2729208 0.1613699 0.47129871 -0.30045520 1.00000000 -0.05151092
## GDD 0.9677481 0.8233183 0.46378408 -0.19039540 -0.05151092 1.00000000Let us keep GDD, bio7 and FDD for plotting NMDS (low autocorrelation, high r2 in environmental fit)
Microclimatic ranges in time and space
…. (Eduardo) (quizás repensar esta parte, pero es importante mantenerla porque es la clave).
DATA AVAILABILITY
The original data, R code for the analysis and creation of the manuscript can be accessed at the GitHub repository https://github.com/efernandezpascual/picos. A version of record of the repository is deposited in Zenodo.
LITERATURE CITED
Figures
Figure 1 Study system.
Figure 1 Study system.
Figure 1 Study system.
Figure 1 Study system.
Figure 1 Study system.
Tables
| Taxon | GDD estimate | GDD p | FDD estimate | FDD p | rho2 | Cold & Frozen | Cold & Snowy | Hot & Frozen | Hot & Snowy | Frequency 2009 | Frequency 2018 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Alchemilla catalaunica | 0.00 | 0.174 | -0.06 | 0.035 | 0.20 | 0 | 17 | 0 | 74 |
|
|
| Androsace villosa | 0.01 | <0.001 | -0.01 | 0.376 | 0.49 | 0 | 0 | 100 | 100 | 17.6 | 17.2 |
| Arabis alpina | -0.01 | 0.004 | -0.01 | 0.396 | 0.22 | 25 | 60 | 0 | 0 |
|
|
| Arenaria grandiflora | 0.00 | 0.004 | -0.01 | 0.272 | 0.18 | 22 | 62 | 0 | 0 |
|
|
| Arenaria moehringioides | -0.01 | 0.001 | 0.01 | 0.03 | 0.34 | 97 | 60 | 0 | 0 | 1.6 | 3 |
| Armeria cantabrica | -0.01 | <0.001 | -0.01 | 0.029 | 0.39 | 65 | 98 | 0 | 0 | 1.6 | 0.2 |
| Carex sempervirens | 0.00 | 0.007 | -0.02 | 0.007 | 0.18 | 1 | 38 | 36 | 97 | 11.2 | 9.8 |
| Erigeron alpinus | 0.00 | 0.972 | -0.06 | 0.035 | 0.19 | 0 | 34 | 0 | 33 | 0.2 | 0.8 |
| Euphrasia salisburgensis | 0.00 | 0.11 | 0.03 | <0.001 | 0.32 | 91 | 1 | 100 | 33 | 4.4 | 3.2 |
| Festuca glacialis | -0.02 | 0.001 | 0.01 | 0.278 | 0.67 | 100 | 100 | 0 | 0 |
|
|
| Festuca hystrix | 0.01 | <0.001 | 0.00 | 0.644 | 0.30 | 0 | 1 | 98 | 99 | 17.1 | 17.9 |
| Galium pyrenaicum | 0.00 | 0.9 | 0.03 | <0.001 | 0.37 | 95 | 2 | 96 | 3 | 0.9 | 0 |
| Helianthemum canum | 0.01 | <0.001 | -0.01 | 0.333 | 0.40 | 1 | 3 | 100 | 100 | 56.2 | 59.2 |
| Iberis carnosa | 0.00 | 0.023 | 0.02 | <0.001 | 0.27 | 98 | 26 | 21 | 0 | 0.2 | 0 |
| Lotus corniculatus | 0.00 | 0.045 | -0.04 | 0.05 | 0.18 | 0 | 10 | 0 | 85 | 0 | 0.1 |
| Scilla verna | 0.00 | 0.006 | -0.05 | 0.035 | 0.25 | 0 | 6 | 0 | 96 | 6.2 | 7 |